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NUMERICAL  SOLUTION  OF  A  SINGULAR 
INTEGRAL  EQUATION  ARISING  FROM  A 
SEQUENTIAL  PROBABILITY  RATIO  TEST 


Sherwood  Sairm 
USAF  Armstrong  Laboratory 


STATEMENT  OF  PROBLEM 


In  a  recently  derived  sequential  probability  ratio  test  (SPRT)  [1]  to  discrimi¬ 
nate  between  the  two  simple  hypotheses:  Ho  ■.  a  =  ao  and  H\  :  cr  =  cri  (0  < 
o'o  <  Cl)  the  relevant  probability  distribution  is  given  by 


9d{z;c)  = 


z  >  -d 


0, 


otherwise. 


where  d  =  In  ^  >  0,  and  c  =  c(cr)  =  y  >  0.  The  operating 

characteristic  (OC),  q,  satisfies  the  integral  equation  [2] 

q{-b-,k,a)=  f  gd{z-,(T)  dz  +  f  gd{z  +  b-,a)q{z-,k,a)  dz,  0  < -b  <  k 

oo  Jo 

which  in  this  case,  with  y  =  —6,  reduces  to 

/•k  1  z+d-v 

q{y,k,(7)  =  fo{y,d,c)  +  Ci  .  - - e~  2c  q{z-,k,(T)  dz,  0  <  y  <  k 

Jko{y,d)  y/ Z  +  a  —  y 

(1) 


where 


koiy,  d)  =  max{0,  y  -  d},  Ci  =  ci(cr)  =  and 


/o(2/;  d,  c) 


2$  -  1,  for  y  <d, 

0,  y>d. 
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Here  $  is  the  distribution  function  of  the  standard  normal  distribution.  The 
integral  is  singular  for  those  values  of  y  satisfying  k  >  y  ^  d. 

For  given  values  of  q:,/9  (0  <  <  l),(ro,  and  cri  (0  <  (Tq  <  ctj),  the 

problem  is  to  find  (unique)  values  (k*,  b*)  of  (A:,  b)  for  which  q(—b*;  k*,  cro)  = 

1  —  a,  and  g(— 6*;  A;*,  (Ti)  = /3. 


METHOD  OF  SOLUTION 


For  given  values  of  (0  <  <  l),o‘o?  ^nd  ci  (0  <  (Tq  <  cr\)i  let 

b  =  bci,ao{k)  denote  the  solution  of  q{—b\  k,  <Jo)  =  1  -  oc,  and  b  =  6/3, (A;) 
the  solution  of  q{-b,k,ai)  =  p.  Then  k*  will  merely  be  the  solution  of 

Once  k*  is  found,  6*  is  immediately  given  by  b*  = 
ba,<roik*).  Thus,  the  basic  problem  is  to  solve  Equation  (1)  for  q{y)  for  each 
given  k  and  cr.  ^ 

By  defining  ^  as  qe~'^,  /o  as  foe~^,  andc2  =  C2(<T)  as  Ci{cr)e  , 
Equation  (1)  becomes 


q{y- k,  a)  =  fo{y,  d,  c)  +  C2  ^ 


ko{y,d)  y/z  +  d-y 


dz,  0  <  y  <  k 


(2) 


With  the  change  of  variable  s  =  k  —  z  m  the  integral.  Equation  (2)  be¬ 
comes 


q{y\  k,  a)  =  foiy;  d,  c)  -|-  C2 

J  0 


k-ko{y,d)  —  5;  cr) 


\/k  -{■  d  —  y  —  s 


ds,  Q  <  y  <_k  (3) 


Or, 


q{y,k,a)  =  fo{y]d,c)  +  C2 


jk  g{k  s;k,<7)  0  <y  <  min{d.  A:} 

-^yk-i-d-y-s 

rk+d-y  g-(k-s;k,cr±  ^  d<y<k 

•^0  y/k+d-y-s  ’  -  ^ 


(4) 
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This  integral  equation  is  solved  numerically.  We  have  used  two  different 
methods.  The  first  method  solves  the  problem  by  interpolating  the  numerator 
in  Equation  (4)  on  each  appropriate  sub-interval  by  a  quadratic  polynomial 
and  then  integrating  the  resulting  integrand  to  obtain  a  set  of  weights  for 
the  numerical  scheme.  Details  of  this  method  can  be  found  in  the  Appendix. 
This  method  is  relatively  easy  to  implement  and  does  not  rely  on  any  special 
software  routines  other  than  a  solver  for  a  system  of  linear  equations. 

The  second  method  makes  indirect  use  of  special  software  routines  avail¬ 
able  commercially.  Here  the  sought  function  q{y)  is  discretized  on  the  inter¬ 
val  [0,fc]  to  yield  the  approximations  of  q{nh),{n  =  0,1,..., A'),  where 
h  is  the  (uniform)  mesh  size  and  k  =  Nh.  Because  the  integral  equation 
(4)  is  only  (weakly)  singular  when  d  <  y  <  k,  the  manner  in  which  it  is 
discretized  depends  on  y.  For  y  <  min{d,  fc},  and  using  Adams  method 
and  Backward  Differentiation  Formulae  method  [3],  one  obtains  a  set  of 
Ni{=  min{nd  =  d//i,  N})  equations  of  the  form 


Pi-1 

qi  =  fo{ih)  +  C2h{ 


WNjqN-j 


N 


UN-jqN-j 


j=o  ”  j)h  j=pi  yj {N  ^  i  -i-  Titi  j)^ 


(5) 


for  i  =  0, 1, . . . ,  Ai.  Here  W  and  w  are  the  appropriate  starting  and  convolu¬ 
tion  weights  of  the  method  respectively  [3],  and  pi  is  the  order  of  the  method 
which  we  took  to  be  4. 

For  y  '>  d  and  using  a  different  Backward  Differentiation  Formulae  method 
[4],  one  obtains  a  set  of  A  —  Ai  equations  of  the  form 


2p2-2  N+Ud-i 

qi  =  foiih)  +  C2\/^{  Y  ^N+nd-i,j^N-j  +  Y  ‘^N+nd-i-j^N-j}  (6) 

j=0  i=2p2-l 

for  z  =  Ai  -f  1, . . . ,  A.  Here  and  are  appropriate  fractional  starting 
and  fractional  convolution  weights  of  the  method  respectively  [4],  and  p2  is 
the  order  of  the  BDF  method  which  we  took  to  be  4.  The  weights  in  both 
Equation  (5)  and  Equation  (6)  are  obtained  from  appropriate  NAG  routines. 
The  two  sets  of  equations  constitute  a  set  of  A  -f  1  linear  equations  in  the 
A-ll  unknowns  solved  using  any  number  of  standard  linear 

equation  solvers.  We  used  one  of  IMSL’s  routines  [5]. 
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RESULTS 


On  the  test  problem  (to  be  described)  as  well  as  the  actual  problem,  the  two 
numerical  methods  have  yielded  similar  results.  Thus  we  will  not  distinguish 
the  two  in  the  following  discussion. 

We  tested  the  numerical  method  above  on  a  simple  problem  whose  ana¬ 
lytical  solution  is  known.  The  problem  consists  of  solving  Equation  1  with  a 
new  /o  given  by 

foiy;  d,  c)  =  e^{l  -  2c2{^Jk  +  d-y-  ^maxlO,  y  -  d}  -  {y  -  d))}  (7) 

The  analytical  solution  for  this  problem  is  simply 

q{y‘i 


or, 


q{y]k,a)  =  1 


The  numerical  result  of  solving  this  simple  problem  is  displayed  in  Figure 

1.  Here  q  is  plotted  against  y.  In  this  particular  test,  the  parameters  are 
A:  =  4,  (To  =  1,0-1  =  2.  q{y,k,a)  was  calculated  twice,  once  for  cr  =  1  and 
once  for  O'  =  2.  Both  of  these  are  shown  in  Figure  1.  As  the  analytical 
solution  of  q  does  not  depend  on  a,  they  are  exactly  the  same,  and  represent 
faithfully  the  anayltyical  solution  q  =  1.  Note,  however,  that  if  we  had  plotted 
q  =  qe'k  —  instead,  the  two  would  definitely  be  different,  since  c  =  c{a) 

is  a  function  of  cr. 

To  illustrate  how  this  method  solves  the  real  problem  (Equation  1),  we 
pick  (Jo  =  1  in  the  null  hypothesis  Hq  and  Ci  =  1.2  in  the  alternative  hypothe¬ 
sis  Hi.  Assuming  a  =  0.05  and  =  0.1,  we  calculated  the  curves  (in  the  k-y 
plane)  qiy;k,(7o)  =  1  -  a  and  q{y;k,(Ti)  =  /?.  These  are  displayed  in  Figure 

2.  These  curves  intersect  at  {k*,y*)  =  (4.75,2.18).  Finally,  in  Figure  3,  we 
plotted  q{y-k*,(To)  and  q{y,k*,ai).  Indeed,  when  y  =  y*,  q{y,k*,(7o)  =  0.95 
and  ^(y.  A:*,  (7i)  =  0.10. 
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CONCLUSION 


In  this  paper,  we  have  proposed  two  numerical  methods  to  solve  a  singular 
integral  equation  which  arose  in  the  study  of  a  certain  Sequential  Probability 
Ratio  Test.  One  method  uses  an  easily-implemented  scheme  and  the  other 
method  involves  casting  the  singular  portion  of  the  problem  into  a  weakly 
singular  equation  of  Abel  type  and  making  use  of  proven  software  for  this 
type  of  equation  ([6]  as  implemented  by  NAG  [4])  to  solve  the  problem.  Both 
methods  have  yielded  similar  results  when  tested  (successfully)  on  a  simple 
problem  and  when  applied  to  the  real  problem. 
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APPENDIX 

A  Simple  Method  for  Solving  Weakly  Singular  Integrals 


In  solving  Equation  (4),  the  main  difficulty  is  solving  a  weakly  singular  inte¬ 
gral  of  the  form 


/(^) 

Ja  (x  —  C)“ 


dx 


(A.l) 


where  a  <  1  and  c  <  a  <  b.  For  the  problem  in  this  paper  o:  =  |.  If  we 

interpolate  f(x)  with  a  (the)  polynomial  P{x)  =  degree  n  at 

the  points  a  =  xq  <  xi  <■■•<  Xn  —  b,  then  the  integral  in  (A.l)  can  be 
approximated  by 


f  ^  ELo 

Ja  {x  —  c)“ 


dx  = 


n  fb 
k=0 


{x  —  c)“ 


dx  =  '^  UkWk  =  w^a, 


k=0 


where 


(A.2) 


Wk 


Ja  (x  —  c)" 


(x  —  cy 

w  =  (wo, . . . ,  u;„)^,  and  a  =  (oq,  . a„)^.  It  is  straight  forward  to  show  that 


Wk  =  Qk(b,c)  -  Qk(a,c), 


where  for  all  y  <  x, 

‘  (') 

Qkix,y)  =  E 

*  +  i  -« 

Since  the  polynomial  P{x)  =  YJk=o  ^kX^  interpolates  f{x)  at  the  n  -fl  points 
a;o, . . .  ,x„,  each  Uk  depends  only  on  and  {/(xi)}”_o.  In  fact,  if  V  is 

the  Vandermonde  matrix  corresponding  to  {x,}"_o,  then 
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Va  =  f, 


(A.3) 


where  f  =  (/(xq),  •  •  • ,  f{xn)f-  Recall  the  Vandermonde  matrix  correspond¬ 
ing  to  {xjf^o  is 

... 

•••  <  / 

Now  by  combining  Equations  (A. 2)  and  (A.3),  an  approximation  to  the  in¬ 
tegral  in  (A.l)  can  be  written  compactly  as 

w^V-^f  (A.4) 

In  solving  Equation  4,  we  have  chosen  the  degree  n  of  the  interpolating 
polynomial  i^(x)  to  be  3.  We  could  have  easily  chosen  it  to  be  any  other 
larger  integer.  By  repeatedly  using  Equation  (A.4)  on  successive  sets  of  n  + 1 
mesh  points  in  a  chosen  partition  {0  =  xo,xi, . . . ,  x;v'  =  of  the  interval 
[0,  k]  (see  Equation  (4)),  we  immediately  obtain  a  system  of  V  -f  1  linear 
equations  in  the  V-fl  unknowns  {/(xo), /(xi), . . .  ,/(xjv)}-  This  system  can 
now  be  solved  by  using  any  standard  system  of  linear  equations  solver. 


/I  Xo  Xq 

1  Xi  Xj 

\  1  Xji 
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Figure  3.  q(y)  at  k*  for  two  sigmas 


11 


